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We investigate the existence and stability of gap vortices and multi-pole gap solitons in a Kagome 
lattice with a defocusing nonlinearity both in a discrete case and in a continuum one with periodic 
external modulation. In particular, predictions are made based on expansion around a simple and 
analytically tractable anti-continuum (zero coupling) limit. These predictions are then confirmed 
for a continuum model of an optically-induced Kagome lattice in a photorefractive crystal obtained 
by a continuous transformation of a honeycomb lattice. 

PACS numbers: 



I. INTRODUCTION 

Recent years have seen a great deal of interest in 
Hamiltonian lattice and quasi-discrete systems due to 
their relevance as models of experiments coming from 
various branches of physics. An early example is that 
of the nonlinear optics of fabricated AlGaAs waveguide 
arrays The interplay of discreteness and nonlinear- 
ity there led to the emergence of numerous phenomena 
that have gathered considerable attention subsequently, 
such as Peierls-Nabarro potential barriers, diffraction and 
diffraction management and gap solitons Js], to name 
just a few. See, for example, the reviews [3, [a] and refer- 
ences therein. 

More recently, there has been growing interest within 
nonlinear optics in the area of optically-induced pho- 
tonic lattices in photorefractive crystals such as stron- 
tium barium niobate (Sra;Bai_a;Nb03, commonly abbre- 
viated SEN) . The original theroetical proposal of [a] was 
followed quickly by experimental realizations [tI, [8|| , and 
the foundation was thus set for the observation of a di- 
verse array of novel nonlinear phenomena in this setting. 
Amon g ot hers, these phenomena include dipole mul- 
tipole [10| , necklac e 111 |, gap and rotary [l3| solitons 
as well as discrete Il3l, Il4| and gap [IB] vortices, higher 
order Bloch modes |lq|. Zener tunneling [13], as well as 
localized modes in honeycomb hexagonal ^9] and 
quasi-crystalline (20| lattices, and Anderson localization 
[2ll | (see, e.g., the reviews [H, for additional exam- 
ples) . 

A considerable effort along these lines has been dedi- 
cated to the recently emer ging area of non-square lattices 

[H, [H, Hi, Hi, in 0, Hilia S [sH, m IM, H [H, Ei]. 

Furthermore, the majority of these studies have dealt 
with the case of a focusing non-linearity rather than a de- 
focusing one. Coherent structures in the latter case have 
received relatively limited attention until very recently, 
for example in the study of fundamental and higher or- 
der gap solitons More complicated gap structures, 
such as multipolcs and complex valued vortices are only 
starting to be explored in square lattices [stI, (sl] . A the- 
oretical framework has been developed in parallel to this 



work, stemming from one-dimensional and square lattices 
[3^, [iO]. However, the predictions of the latter can also 
be translated to contours (or paths) in non-square ge- 
ometries [2^ [3l| , based on arguments of dimensionality 
reduction along the contour. 

In this work we will focus on the so-called Kagome lat- 
tice, which is encountered often in nature and has a very 
rich structure. In the solid-state community and other 
areas of physics and science these lattices and many oth- 
ers have been explored for decades [4l[ , but are becoming 
more prevalent recently (4^ . l43j . Furthermore, low tem- 
perature properties of atomic quantum (ultracold Bose 
and Fermi) gases have been studied in the trimerized 
Kagome lattice [4l| . 

Motivated by recent advances in optically-induced lat- 
tices in SBN, we will explore a Non-linear Schrodinger 
(NLS) model in both its discrete (DNLS) manifestation 
as a set of difference equations adhering to the symmetry 
of the lattice and modeling coupled oscillators, and in the 
analogous continuum setting using a partial differential 
equation with an external potential having the appropri- 
ate symmetry. In particular, since the continuum model 
is motivated by experiments with SBN, the nonlinearity 
will be saturable [13, H^. We will investigate prototypical 
contours (or paths) of localized structures in this lattice, 
consisting of six sites as well as four sites, and being both 
real, and complex valued with continuous phase (modulo 
2tt). 

Our main findings in what follows are that 

• Certain structures are stable, such as the in-phase 
gap hexapole and single-charge six-site gap vortex 
on the honeycomb cell, and the in-phase/out-of- 
phase quadrupole on the "hourglass cell" (see Fig. 
[1]). Other configurations are unstable. 

• In the continuum model, continuations of solutions 
in the first band-gap pass through the second band 
as quasi-localized structures and then become fully 
extended in the second band-gap. However, discon- 
tinuous extensions, i.e. new continuations of the lo- 
calized structures, are found to exist in the second 
band-gap simultaneously with the extended states. 
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• The result of the evolution of the dynamical in- 
stability in these lattices is more complex than in 
the square lattice case, and may involve not only 
degeneration to single-site solitons but possibly to 
multi-site solitary wave structures, and, in the dis- 
crete case, often the formation of robust breathing 
states, consisting of multiple sites (possibly even as 
many as in the original configuration). In fact, we 
have found some clear breather formations recur- 
ring in multiple simulations: 

— Two nearest-neighbor or opposite sites in- 
phase with each other and with oscillating am- 
plitudes of comparable magnitude. 

— Two next-nearest-neighbor sites out-of-phase 
with each other and with oscillating ampli- 
tudes of comparable magnitude. 

— Two nearest-neighbor sites having different 
amplitudes and oscillating between the same 
phases and opposite phases depending on 
whether the amplitudes are further from or 
closer to each other, respectively. 

And, in the continuum version, either all or most 
of the initially populated wells remain populated 
for long propagation distances, with the instability 
manifesting itself only as phase reshaping. 

The presentation will be structured as follows. In sec- 
tion |TT] we provide the setup of the problem, the back- 
ground and the theory. Then in section IIIII we will sys- 
tematically explore the relevant numerical results. Fi- 
nally, in section IIVI we will summarize our findings and 
present our conclusions. 



II. SETUP 



where ej is a translation by one site in the positive di- 
rection along J and J is one of a site-dependent subset 
of two of the principal lattice vectors = (l,\/3)/2, 
be; — (— 1, V3)/2, or Cd = (1,0) of the discrete Kagome 
lattice presented in Fig. [1] and e is the coupling between 
sites. The non-linear term is taken to be a cubic Kerr 
non-linearity as follows 



U{^,\U\')^-\U\' 



(3) 



The simulations for the static results in the discrete 
model were performed in the domain Dh\K, where 
Dh = [1, . . . , 33] X [1, . . . , 33] is the discrete lattice do- 
main corresponding to a triangular lattice and K = 
{{2m + l,2n+ l)\{m,n) € [0,16]^}. For dynamical evo- 
lution, the solutions were buffered with 40 (or more) 
nodes on all sides to prevent radiation scattering from 
the boundaries. 




The description of the setup of the problem, includ- 
ing the background and theory is organized into three 
sections: first, the preliminary material III Al second the 
existence considerations ITl Bl and third the stability con- 
siderations III CI 



FIG. 1: (Color online) The discrete Kagome lattice structure 
is presented above. The six-site contour is given by the blue 
circles, while the four-site "hourglass" contour is given by the 
green squares. 



Preliminaries 



The continuum version consists of defining C 
where is the two-dimensional Laplacian on M 



R and 



We introduce the following complex-valued non-linear 
evolution equation 



AA(x,|t/p) = 



l+/(x) + |C/|2' 



(4) 



-lU, = [C+J\f{^,\U\^)]U, 



(1) 



where U is a function of z e IR.+ and x G Z x Z in 
the discrete version or else x e R x M in the continuum 
version. First, we will consider the discrete version with 



C 




(2) 



with 



/(x) = lo \fiiii)e'''^'-'' + e"'''"^-^ + e^'^'^^-'^l (5) 

the optical lattice intensity function formed by 
three (p=0) or four laser beams with /i(x) = 
gifepx/(i+4p/3) cos[pkx/{l + 4:p/3)],hi = (l/(l + 4p/3),0). 



- (^2(l+4p/3) ' ~"^)' ^3 - (^2(l+4p/3) ' 
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p — > 3/2 the lattice transforms from the well-known hon- 
eycomb interference pattern into the richer Kagome lat- 
tice. The latter lattice features both the hexagons from 
the honeycomb lattice and the equilateral triangles from 
the triangular one, and each node has four neighbors sim- 
ilar to the square lattice. 

Here /q is the lattice peak intensity, z is the propa- 
gation constant and x = ix^y) are transverse distances 
(normalized to = 1 mm and Xs — Vs — l/^ui), 
is proportional to the applied DC field voltage, D — 
z^Xj {^'Krif.Xsys) is the diffraction coefficient, A is the 
wavelength of the laser in a vacuum, d — ^-n jk is the pe- 
riodicity in the x-direction [dj 1/3 is the periodicity in the 
y-direction) and Ue is the refractive index along the ex- 
traordinary axis. We choose the lattice intensity /q = 1, 
and (d, i?o, A, ne)=(90,8,532 nm,2.35), consistent with a 
typical experimentally accessible situation [H, [sj- A 
plot of the potential intensity field created by the opti- 
cal lattice is shown in Fig. [5] to illustrate the locations 
where our localized configurations will live. The non- 
dimensional value D — 18.01, and we note that this dis- 
persion coefficient is equivalent to rescaling space by a 
factor \fD as, e.g. in [4y|. 

The numerical simulations are performed in a rectan- 
gular 120 X 120 grid corresponding to the domain size 
4(i X 8(i/\/3 (i.e. four periods of the lattice in each di- 
rection), using a rectangular spatial mesh with Aa: = 1.5 
and Ay « 1.732. Regarding the typical dynamics of a so- 
lution when it is unstable, we simulate the z-dependent 
evolution using a Runge-Kutta fourth-order scheme with 
a step size Az = 0.01. 



B. Existence considerations 

Assuming a stationary state u exists, and letting the 
propagation constant — represent the (nonlinear) real 
eigenvalue of the operator of the right-hand-side of Eq. 
((!]), the corresponding eigenvector u is a fixed point of 
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FIG. 2: (Color online) A spatial (x-y) contour plot of 
the effective potential created by the ordinary polariza- 
tion standing wave [lattice beam in Eq. ((5]l]. Points 
A, B, C, D, E, F, G and H define the relevant potential min- 
ima for the various configurations we will consider. The 
contour {A, B ,C, D, E, F} is the honeycomb cell, which can 
be considered to tile part of the lattice. The set of sites 
{B, F, G, H} comprise the "hourglass" cell contour we will 
consider. Together with A, these sites comprise another cell 
which tiles the remaining part of the lattice. 
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[H + C+U{x,\uf)]u = 0. 



(6) 
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In the discrete case, we performed a continuation in 
e. In particular, if one indexes the sites by (m,n), then 
solutions in the limit e — *■ can easily be found of the 
general form Um.n = y^exp {— i(/it — dm,n)} for any ar- 
bitrary 0„i,n G [0, 27r) We can linearize Eq. ^ 
around the solution for e = denoted by uq, accounting 
for complex valued perturbations by considering the con- 
jugate of Eq 
Jacobian of 1 
of £, is 



3]) as well, which has the solution Uq. The 
for e = 0, or equivalently, in the absence 



FIG. 3: (Color online) The discrete in-phase hexapole solu- 
tions are presented. In the first two columns the profiles (left) 
and linearization spectra are given before (top, e = 0.061) 
and after (e = 0.085) the first Hamiltonian Hopf (HH) bifur- 
cation. The top right panel depicts the theoretical predic- 
tions of the linearization eigenvalues bifurcating from the AC 
limit (dashed) as well as the actual numerically computed 
ones (solid). The bottom right panel is (see Eq. (|10|l '). 
shown on a log scale, where we can observe the decrease in 
the effective power, as the coupling strength increases. 



J{u) = [fi + d{J\fu, [Afu]*)/d{u, u*)]. 



(7) 



We may also take the coupling e, when sufficiently 
small, as the small parameter in the expansion with 



[u,u*] = [uo,Uq] + sui. If we denote by the op- 
erator C for coupling e, then the first order correction in 
e to Eq. ® is 
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FIG. 4: (Color online) The dynamics of the solution given in 
the bottom row of Fig. |3] is presented. The top row shows 
snapshots of the modulus for various z, while the next row 
shows the individual amplitudes at the relevant sites. The 
structure survives for a while but ultimately disintegrates, 
due to the instability, into two populated nearest-neighboring 
sites whose amplitudes breathe closer and further from one 
another while the phases oscillate between opposite and same, 
respectively (see the bottom two rows). 
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Projecting this map onto the kernel of J^{uq) elimi- 
nates the first term and we are left with the condition 



where we use (•,•) the standard inner product on the 
Hilbert space P. We let /i = 1 without loss of gener- 
ality and denote by j the indices (m, n) along the one- 
dimensional contour. The non-trivial part of the Jaco- 
bian J{uq) decouples into a direct sum of 2 x 2 blocks if 
there are N excited sites in the contour. For each j there 
is a nontrivial element (e*^^ , — e"*^^ )^ € ker{ J, (mo)}. So, 
the condition for existence of solutions to Eq. ^ with 
e > reduces to the vanishing of the vector function g(0) 
of the phase vector 6 = {9i, . . . ,6m) where 



FIG. 5: (Color online) The continuum in-phase hexapole is 
presented in these panels. The top two left panels show the 
power, P (top) and instability growth rate (bottom) as given 
by the maximum real part of the linearization spectrum. The 
first band is given to the right, beyond which is the semi- 
infinite gap (it is displayed wider than it actually is for visi- 
bility, because its actual width is narrower than a pixel at this 
scale) , the second band is in the middle, and the third band is 
at the far left. The blue branches in the second gap are actu- 
ally discontinuous extensions of the localized modes from the 
first gap, which collide in a saddle-node bifurcation and disap- 
pear as can be observed in the inset panels in the upper right 
corners (this is consistent throughout the following images, 
but the closeups will not be shown). Solutions marked on 
these plots with the letters a,b,c,d,e,f,g and h are presented 
in the remaining panels. The second and third columns of 
the top set display the principal two solutions a and b, re- 
spectively, with full panels of their linearization spectra be- 
low them. The bottom six panels have miniature sub-panels 
with the corresponding spectra embedded in them. There are 
stable first (a), and also second (g and h) gap soliton struc- 
tures. The solitons (c,d), with energy in the second band, are 
unstable. 
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FIG. 6: (Color online) The evolution of the unstable solu- 
tion given in Fig. [S] (b) is shown. The phase is shown as 
arg(u)x{(2;, > O.Sniaxfj;^^) at various times be- 

cause the original configuration is preserved for a very long 
propagation distance (x is the indicator function which anni- 
hilates the field outside of the set on which it is defined) . The 
relative phases of the configuration break up after z — 100. 



(8) 



subject to periodic boundary conditions. Wc consider 
primarily contours M within the subcategory of those 
for which \0j+i — Oj\ — A0 is constant for all j € M, 
|6»i - ^iMil = AO and Ae\M\ = mod 27r. A standard 
Newton fixed-point solver is used to construct branches 
of solutions to Eq. ([5]) in e from the AC limit. 

For the continuum problem where x g M^, there ex- 
ists no such analytical solution from which to construct 
a continuation. On the other hand, it is well-known that 
localized solutions exist for values of the propagation con- 
stant /i in the complement of the linear spectrum (ie. the 
so-called spectral gap) defined by the following eigenvalue 
problem (also known as the linearization around the zero 
solution), 



[/!-£- 7V(x, 0)]u = 0. 



(9) 



These solutions are exponentially localized in space, 
so-called gap-soliton, states of the original nonlinear par- 
tial differential equation. Since the parameter of inter- 
est is /i, diagnostics are plotted against /i. Continua- 
tions in this parameter can be found with a fixed-point 
solver and an initial guess of a collection of Gaussian 
wave-packets in the appropriate configuration. Using a 
standard eigenvalue solver package implemented through 
MATLAB, we identified the first two spectral gaps for our 
given parameters and grid size to be Gi ~ (3.545, 4.9454) 
and G2 ~ (2.178,3.515). It is worth noting that the 
bands and gaps remain very close to the same widths 
for much smaller discretizations (ie. much larger grids). 
For instance, with 300 nodes in each direction we have 
G2 = (2.125, 3.463), so the change in width of the band- 
gaps is an order of magnitude closer to convergence than 
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FIG. 7: (Color online) The same panels as Fig. [3] except for 
the unstable out-of-phase hexapole. The top row solution is 
for e — 0.05, while the bottom one is for e = 0.16. 



the position (modulo translation). We use the bands ap- 
propriate to the discretization in order to compare them 
with the bifurcation structure of solutions. 

The localized states u of the continuum version of ^ 
were obtained using the Newton-Krylov fixed point solver 
nsoli from 47], which utilizes a GMRES iterative algo- 
rithm, based on residual reduction in successive Krylov 
subspaces, in order to minimize the memory necessary 
for the linear solver within each step of the Newton al- 
gorithm. Some care has to be taken to handle the large 
size of the representation of a 2D continuum domain. 
A pseudo-arclength continuation 48] was used to follow 
each branch and locate the bifurcations which occur at 
the edges of the bands. 

The square root of the optical power of the localized 
waves is defined as follows: 



P 



\U\^dS 



1/2 



(10) 



where in the continuum problem, dS = dxdy, while in the 
discrete problem we define the corrsponding sum (divided 

by V^). 



C. stability considerations 

Stability is examined by linearizing Eq. ([T]) and its 
conjugate around an exact stationary solution (u,m*)^ 
to Eq. If we assume the perturbation is separable 

of the form u = e^'^w{x), then we have the following 



eigenvalue problem for i\ 



[J{U)+C2\ 



> V 



0, 



(11) 



where J is defined in Eq. ([7]) and £2 , 12 are the 2x2 block 
diagonal matrices with C, I (respectively) on the diago- 
nal. The eigenvalue problem imposes no restriction on 
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FIG. 8: (Color online) The same as Fig. |4]for the out-of-phase 
solution given in the top row of Fig. [7] The original configura- 
tion clearly shifts very rapidly, although the different sites re- 
main populated for a long propagation distance. The largest 
amplitude pair which is separated by one node (i.e. next- 
nearest neighbors) remains very close to exactly out-of-phase, 
while the next smallest, which is also next-nearest neighbor 
to both, passes from the phase of one to the other. The other 
three sites do not exhibit any phase correlation, although at 
times two have matching phase and are out-of-phase with the 
other. 



the eigenvectors v. This form of the Jacobian is identical 
for the discrete and continuum versions of the problem, 
up to the definitions of the operators A/" and £ and the 
domain in which the spatial variables live. The modest 
size of this matrix for the discrete case, [2 x 33 x 33]^ 
is not a problem for a full diagonalization and we im- 
plement the MATLAB function eig to do so. On the 
other hand, the matrix for the model of the continuous 
domain is [2 x 120 x 120]^ and cannot be inverted. For- 
tunately a standard finite difference discretization leads 
to a sparse banded matrix which is perfectly suited for 
Arnoldi iterative algorithms which minimize memory and 
use successive approximations of eigenvalues and eigen- 
vectors until convergence. Such a method is implemented 
by the MATLAB function eigs, which we use here. 

Twofold symmetry over each of the real and imaginary 
axes is guaranteed by the fact that the Jacobian of the 
full problem, which defines the linearization at any point, 
H has the property that JHJ = (where J is the 
canonical symplectic matrix having the properties J J = 
— / and JJ^ = /, which implies = M is symplectic or 




FIG. 9: (Color online) The same panels as Fig. (5] but for the 
out-of-phase hexapole. 
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FIG. 10: (Color online) The same panels as Fig. |3] except for 
the stable single charge vortex solution and the modulus of 
the profiles are given, i.e. |up instead of u. The particular 
solutions are for e = 0.06 (top) and e = 0.11 (bottom). 



7 




-0. 01 0.01 



arg(u) 



arg(u) 



-inn 












o(u),a 



a(u),b 



-2 

-0.01 



arg(L) 



0.005 
Re(?.) 



-0.01 0.01 





Re(?.) 



|u|^,e 



arg(L) 



{j7(uo)}, which are both mapped to zero eigenvalues of 







{i7(mo)}- So, one can foUow the same proce- 



dure for the eigenvalue problem of the linearization (Eq. 
Ilip of the original problem, except projecting also onto 
the generalized kernel and expanding to lower order {^/e) 
and, hence, obtain a mapping between its eigenvalues and 
those of the Jacobian of g. In particular, for each eigen- 
value 7j of A4 , the full linearization around a stationary 
solution with non-zero nodes in M, given by Eq. (jlip . 
will have an eigenvalue pair Aj given, to leading order, by 
Xj — ±^27j£ in the case that the sites in M are nearest 
neighbors of each other. If the excited nodes are sepa- 
rated by a site then, for the DNLS model eigenvalues, e 
is replaced by in the previous relation. The Jacobian 
matrix A4 has the following form: 




k, 

k±l, 



(12) 



- 0.100.1 

in^ 



i| If I • ^ 



|u|^,d 



arg(L) 



Re(i) 
-0.1 0.1 



|u|=,f 



arg(L) 



-100 



i 

o 1 ° 

i 



















-50 



50 



-50 







50 



FIG. 11: (Color online) The same panels as Fig. O but for the 
single charge vortex solution and again, as in the analogous 
discrete case, the modulus is given in lieu of the field itself. 
Here there are small embedded panels in the top right corners 
of the profile images with the phase of the solution. Also, we 
were not able to identify continuations of the branches from 
the first gap into extended solutions in the second gap in this 
case, although we did identify semi-extended states. 



M^JM — J), and so Re(Aj) ^ implies an instability. 

In the discrete case, the stability will be compared to 
analytical results for small e based on Lyapunov-Schmidt 
analysis of the expansion of the equation around the 
AC limit (see, e.g. and references therein for de- 

tails). For the contour M, there are \M\ eigenvalues 
7j of the \M\ X \M\ Jacobian Mjk = dgj/dOk of the 
diffeomorphism given in Eq. ([8|). Now, for each ex- 
cited site, there are a zero and a negative eigenvalue of 



We now briefly discuss the princiml stability conclu- 
sions, for the defocusing case of [1^, [3§|, which we expect 
to remain valid in the present configuration. Nearest 
neighbor excitations in the defocusing case correspond 
to nearest neighbor excitations in the focusing case, but 
with an additional tt phase in the relative phase of the 
sites added by the so-called staggering transformation 
[39| . Therefore, the in-phase nearest neighbor configura- 
tion in the defocusing case corresponds to an out-of-phase 
such configuration in the focusing case (and should thus 
be stable) On the other hand, next nearest neighbor 
out-of-phase defocusing configurations would correspond 
to next nearest neighbor out-of-phase focusing configura- 
tions and should also be stable (at least in some param- 
eter regimes). By the same token, out-of-phase nearest 
neighbor, and in-phase next nearest neighbor structures 
should be unstable. Finally, the single-charge vortex and 
in-phase hexapoles should be stable, while the double 
charge vortex and out-of-phase hexapoles should be un- 
stable. However, notice that, as discussed in [s^, the 
multipole structures characterized as potentially stable 
above will, in fact, typically possess imaginary eigenval- 
ues of negative Krein signature (see e.g. [49| and ref- 
erences therein). These may lead to oscillatory insta- 
bilities through complex quartets of eigenvalues, which 
arise by means of Hamiltonian-Hopf (HH) bifurcations 
(50| emerging from collisions with eigenvalues of oppo- 
site (i.e., positive) Krein signature. These conclusions 
will be discussed in connection with our detailed numer- 
ical results in what follows. 
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III. NUMERICAL RESULTS 

Now, the above theoretical predictions will be matched 
against systematic numerical simulations. First, for the 
discrete model, we will perform continuations in the cou- 
pling parameter from the AC limit in order to compare 
the resulting relevant eigenvalues from the linearization 
spectrum with the corresponding prediction. Next, we 
will test these results against the continuum model. The 
stability results from the discrete case are expected to 
hold in the sense that there will either be real eigenvalue 
pairs in the spectrum of the solutions whose discrete ana- 
log is unstable close to the AC limit, or else there will be 
intervals of stability and quartets of eigenvalues due to 
HH and inverse HH bifurcations. The continuum model 
reveals not only gap soliton solutions in both the first 
and the second gaps, but also solitons from where the 
branch of solutions from the first gap passes through the 
first band and subsequently becomes extended. Unsta- 
ble solutions were evolved in time in order to observe 
their dynamical behavior. Most solutions in the discrete 
case decompose into breathing configurations with fewer 
populated sites and some interesting phase correlations. 
In the continuum case, most configurations survive for a 
long propagation distance, with instabilities manifested 
only as phase reshaping. 

This section will be composed of two parts, the first 
of which will address configurations with six neighbors 
on the hexagonal cell, the results of which are consistent 
with recent results in the continuum honeycomb defo- 
cusing case [1^ and both honeycomb and hexagonal [Slj 
focusing cases (translated with the appropriate stagger- 
ing transform along the contour). In the second part we 
will look at quadrupoles along the four corners of the 
"hourglass" cell which is unique to the Kagomc lattice. 

A. Vortices and hexapoles in the hexagonal cell 

The results for the six-site configurations in the hexag- 
onal cell are presented in this section. First we will con- 
sider the real-valued configurations of A6' = and — 
TT and then the complex- valued ones where = 7r/3 
and M = 2TT/i. 

1 . In-phase 

Here we will consider the results of the predictions from 
the previous section for the six-site in-phase configuration 
(i.e. /S.9 = 0). This configuration has been predicted to 
be stable. In the discrete model close to the AC limit, Eq. 
[12] predicts, to first order in e, two double pairs of eigen- 
values i^/2e, i\f%e and single pairs at i\f%e and 0. Here we 
digress slightly to discuss the bound of the phonon band. 
We consider plane waves of the form w ~ g^ip^+i^^) ^ in 
each of 2 (of the 3) principal directions, which are used to 
index the 2-dimensional lattice [2^, (m,n). Since there 



are 3 types of nodes in this case, each having neighbors in 
2 of the 3 principal directions (that are the same for the 
hexagoinal lattice), we must consider a linear combina- 
tion of equal 1/3 weights of the corresponding dispersion 
relation for each type. This is equivalent to 2/3 of the 
dispersion relation of the hexagonal lattice, i.e. 



C^w = (4 — -[cos(g) -I- cos(p) + cos(p + q)])e < 6e. (13) 



So, the smallest eigenvalue of the phonon band is given 
by i(l — 6s), and upon its collision with the eigenvalues 
which bifurcated from the origin, a cascade of HH bi- 
furcations ensues. The numerical results are presented 
in Fig. [3l The left column displays two solutions, be- 
fore (top) and after (bottom), the continuous spectrum 
intersects with the bifurcation eigenvalues. The middle 
column has the corresponding linearization spectra. The 
top right panel depicts the imaginary component of the 
bifurcation eigenvalues as given numerically (solid line) 
and by the first order theoretical approximation (dashed 
line). 

The dynamical evolution of the unstable solution from 
the bottom row is displayed in Fig. (¥] Eventually the in- 
stability manifests itself and the original configuration is 
destroyed. Two nearest-neighbor sites remain with dif- 
ferent amplitudes that oscillate. When they are closer 
in amplitude they are out-of-phase while when they are 
further apart, they are in-phase. It is worth noting here 
that a similar phenomenon was found in the hexagonal 
as well as the honeycomb lattices with a focusing nonlin- 
earity in [3l|, except with the relative phases reversed, 
i.e. the sites were in-phase when closer and out-of-phase 
when further apart, presumably due to the nature of the 
nonlinearity (focusing versus the defocusing one here). 

Next we investigate the in-phase hexapole in the con- 
tinuum setting. The solution is stable in the entire first 
band-gap (see Fig. O top two rows and a). When it 
reaches the first band it collides with an unstable branch 
(b) which has two neighboring wells populated out-of- 
phase and disappears in a saddle- node bifurcation. When 
these branches reach the second band, they immediately 
become unstable as they reshape into extended solutions 
(see e,f). There does exist a second-gap soliton solution 
which is stable, (g). This solution disappears in a saddle- 
node bifurcation with a marginally stable solution that 
has next-nearest-neighbor wells on four sides populated 
with intra-site dipoles. The evolution of solution (b) is 
represented in Fig. [6]by the phase of the sites via the fol- 
lowing quantity aTg{u)x{{x,y)\\u\'^ > 0.5max(3,_y) [jup]}, 
where x is the indicator function of the set that annihi- 
lates the field outside that set. All sites remain for a long 
propagation distance, although the phase correlation is 
lost hy z = 300. 
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2. Out-of-phase 

In this section, we consider the unstable out-of-phase 
hexapole {A9 — n). First, in the discrete case presented 
in Fig. [7] the theoretical prediction of Hnearization eigen- 
values, which are exactly a factor —i times those for the 
in-phasc solution, are confirmed for small e. The dashed 
line 1 — 6e, which represents the smallest phonon eigen- 
value, is included here to show that the actual lineariza- 
tion eigenvalues remain bounded by this line, similarly 
to what was observed for square lattices in [39|. A so- 
lution is shown for small coupling and large gap, as well 
as one when the gap is closing and the solution decay- 
ing. The dynamics in Fig. [5] reveal that all sites sur- 
vive for a long propagation distance, and, while there 
is no clear phase correlation between all sites, there is 
between some. For instance, the largest amplitude two 
next-nearest-neighbor lobes remain out-of-phase. This 
is again consistent with a feature that was recently ob- 
served in [sij . since next-nearest- neighbor interactions 
are expected to be the same for focusing and defocusing 
non-linearities due to the staggering transformation [40| . 
The middle amplitude site next-nearest to both of these 
oscillates between in-phase with one and then the other, 
while there is no apparent correlation of the other three 
smaller amplitude sites. 

The same panels as the in-phase case. Fig. [Sj are shown 
for the continuum version of the out-of-phase hexapole 
in Fig. O The solution this time actually collides with 
a four-well structure, which is slightly more stable, due 
to fewer unstable pairs of populated wells. Again there 
are continuations through the second band to extended 
states and again there exists a disjoint branch of sec- 
ond gap states. All states here are unstable. Under dy- 
namical evolution, again all six sites remain for a long 
propagation distance (not shown). However, the phase 
correlation breaks down as early as z = 20, due to the 
instability. 



3. Single charge vortices 

Next we look at the stable single charge vortex solu- 
tion {A9 = 7r/3). The discrete problem is predicted to 
be stable with double pairs of eigenvalues at ±iy/e and 
±i^/3e, and single pairs at ±2iy/e and 0. The prediction 
is confirmed in Fig. [TOland there is good agreement until 
the HH bifurcations set in with the continuous spectrum 
when y/e ~ Ai = 1 — Ge. A cascade of such bifurca- 
tions follows and an example profile and spectrum after 
this time are shown in the bottom left. The evolution of 
the unstable solution from the bottom row was investi- 
gated (not shown) and four sites decompose into essen- 
tially background radiation, while two cells opposite to 
one another inherit most of the power and remain close in 
amplitude and in-phase for a long distance. This is again 
in agreement with the results of [31| , since in-phase oppo- 
site sites on the honeycomb cell are next-to-next-nearest- 
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FIG. 12: (Color online) The same panels in the left two 
columns as Fig. |3]except for the unstable in-phase quadrupole 
in the hourglass cell. The particular solutions shown are for 
£ = 0.03 (top) and e = 0.067 (bottom). 



neighbors, so with a defocusing nonlinearity it is equiva- 
lent to an out-of-phase pair with a focusing nonlinearity. 
Many comparable amplitude out-of-phase breathers were 
found in [3]| . 

The continuum version of this configuration given in 
Fig. [11] (a) collides in a saddle-node bifurcation with an 
unstable configuration with additional populated sites on 
the perimeter, (b). The second gap version (e) collides 
with an unstable state (f) having two intrasite dipoles 
populated outside the original vortex. However, this 
state appears to stabilize closer to the third band. In the 
dynamical evolution of (b, not shown) again the origi- 
nal configuration of "mass" survives for a long propaga- 
tion distance, while the phase correlation decomposes by 
z = 100. 



4- Double-charge vortices 

The last of the six-site configurations we consider is 
the double-charge vortex {A9 — 27r/3). The stability 
predictions are again confirmed for small e; however, due 
to the instability of the branch throughout its existence 
range, we do not present the details here for the sake 
of brevity. The continuum model admits similar branch 
structure for this solution as for previous cases. Here 
also, the solution is unstable for all the cases examined. 
Its dynamics showed the principal sites surviving for long 
propagation distances, although the structure eventually 
disintigates. Again due to the generic instability of this 
branch, numerical details are omitted here. 



B. Quadrupoles in the hourglass cell 

Next, we consider configurations around the outer 4- 
site contour of the hourglass cell. Now, in this case, 
the five sites comprising the hourglass are not a simple 
curve, i.e. the curve crosses itself, and since some nodes 
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FIG. 13: (Color online) The same as Fig. |4]for the in-phase 
quadrupole given in the bottom row of Fig. 1121 The center 
site becomes populated around z — 100 and the remaining 
in-phase tripole (expected to be stable) persists for a long dis- 
tance with two sites almost exactly in-phase, and with equal 
amplitudes, while the other has much larger oscillations in 
amplitude opposite to the other two, but has very similar 
phase. 
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are nearest-neighbors, while others are next-nearest, the 
analysis must be taken to second order for accurate pre- 
dictions of all the bifurcating eigenvalues. Instead, we 
consider only nearest neighbors analytically and extend 
previous results about higher-order interactions to make 
qualitative predictions there. 



C. In-phase quadrupoles 

First, we consider the in-phase quadrupole. The pre- 
diction of the discrete model to first order, i.e. for near- 
est neighbors, is that this configuration will have double 
eigenvalue pairs at zL2iy/e (due to the in-phase nearest 
neighbors predicted to be stable). On the other hand, 
next-nearest neighbors which are in-phase are expected 
to be unstable [23, and indeed a real pair does bifur- 
cate as well. The real pair comes at higher order because 
it is a higher-order splitting. The results are presented 
in Fig. [T2I The panels are the same as for Fig. [3l The 
evolution in Fig. 1 131 reveals a unique structure with three 
sites very close in phase and intensity, two nearly identi- 
cal and one slightly different with its intensity oscillating 



FIG. 14: (Color online) The same panels as Fig. [5] but for 
the in-phase quadrupole. A more unstable configuration with 
the center well populated out-of-phase collides with this one 
close to the first band-edge. 



z=10 z=60 z=80 




FIG. 15: (Color online) The same as Fig. |6] but for the in- 
phase quadrupole from Fig. [14] (b) . The initial relative phase 
loses the correlation after z — 60 and the structure again 
persists for a long distance. 
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FIG. 16: (Color online) The same panels as Fig. [12] except 
for the unstable out-of-phase quadrupole. The particular so- 
lutions shown are for e = 0.03 (top) and e = 0.16 (bottom). 
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FIG. 18: (Color online) The same as Fig. |S] but for the un- 
stable out-of-phase quadrupole given in Fig. [17] (a) . 
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FIG. 19: (Color online) The same panels as Fig. |3] except for 
the stable in-phase/out-of-phase quadrupole. The particular 
solutions shown are for e = 0.03 (top) and e = 0.144 (bottom). 
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FIG. 17: (Color online) The same panels as Fig. [S] but for 
the out-of-phase quadrupole. It collides with a branch that 
has a similar phase pattern as the original configuration. 



FIG. 20: (Color online) The same panels as Fig. [S] but for 
the in-phase/out-of-phase quadrupole. 
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with larger amplitude opposite to the others. One of the 
sites is the originally unpopulated center site, which be- 
comes populated when the other two disintegrate around 
z = 100. 

The continuum version of this configuration disappears 
at the first band edge when it collides with a more un- 
stable solution having the center site populated out-of- 
phase to the others (see Fig. The main branch is 

only weakly unstable from the higher order interactions 
and, in fact, the second-band solution actually becomes 
stable for /i ^ 3. In the dynamical evolution of the struc- 
ture from Fig. [T3] (b), the unstable five site in-phase 
structure remains very robust for a long distance with 
reshaped phase; see Fig. [131 

D. Out-of-phase quadrupoles 

Next we consider the out-of-phase quadrupolc. The 
predicted eigenvalues are the same as those of the in- 
phase case, multiplied by i. They are fairly accurate for 
small £ as one can see in the top right panel of Fig. [111 
The dynamical evolution of the solution given in the top 
row of Fig. [121 reveals two pairs of uneven amplitude 
breathers with phases and amplitudes oscillating oppo- 
site to each other as in Fig. [31 (not shown). 

This solution in the continuum version, as seen in Fig. 
[T71 is always unstable. It collides with a structure that 
has a similar phase pattern, but which is surrounding 
rather than including the original configuration. At the 
point of bifurcation the common structure they share is 
two rows of opposite phase. Again the unstable configu- 
ration persists for a long propagation distance, suffering 
merely a reshaping of the relative phase (see Fig. [T51) . 

E. In-phase/Out-of-phase quadrupoles 

Finally, we turn to the quadrupole solution which has 
its nearest-neighbors in-phase and the next-nearest ones 
out-of-phase. The theoretical prediction for the discrete 
model based on the set of all possible dipole configura- 
tions always implies stability. Indeed the precise first 
order calculation predicts this configuration will also be 
stable with two pairs of eigenvalues at ±12-^. Moreover, 
previous results pol . [39| predict that next-nearest neigh- 
bors which are out-of-phase will be stable, and all those 
in this configuration (except the ones which are also near- 
est) are out-of-phase. The agreement is very good again 
as given in the top right panel of Fig. [TH 

The dynamical evolution of this configuration (not 
shown) again reveals the usual in-phase to out-of-phase 
uneven intensity breather pair, as shown first in Fig. [H 

The continuum version is presented in Fig. 1201 Stable 
first and second band versions of the solution are identi- 



fied and again there are bifurcations at the first and sec- 
ond bands, and also intermediate as well as extended so- 
lutions. The solution that collides with the main branch 
at the first band-edge (b) was propagated (not shown) 
and again the original sites persist and the relative phase 
reshapes after z = 50. 



IV. CONCLUSIONS 

In conclusion, we have presented results of prototypi- 
cal contours (i.e. hexagonal and "hourglass" ) of localized 
structures in a Kagome lattice symmetry with a defocus- 
ing nonlinearity for both a discrete model and an anal- 
ogous continuum model that is relevant to experiments 
on the nonlinear optics of photonic lattices in photore- 
fractive crystals. We have identified stable configurations 
such as the in-phase hexapole and single-charge six site 
vortex on the honeycomb cell, as well as the four-site in- 
phase/out-of-phase quadrupole and the second-gap in- 
phase quadrupole on the hourglass cell. Many of the 
structures admitted not only second gap localized struc- 
tures but also semi-localized structures having energy 
within the second band and extended structures in the 
second band-gap. 

The unstable solutions were evolved in time and all of 
the discrete structures persisted, with recurring dynami- 
cal breathing configurations reappearing in several cases, 
such as the comparable intensity in-phase pair and the 
uneven oscillating intensity pair which are out-of-phase 
when their intensities are closer and in-phase when they 
are further. In the continuum model, there was very 
little deviation from the amplitudes of the initial con- 
figurations and instability was manifested through phase 
reshaping. This suggests that the Kagome lattice is a 
robust structure in which to perform experiments. 

This work suggests many interesting directions to pur- 
sue, the most obvious of which is the experimental real- 
ization of these structures in optical lattices in photore- 
fractive media and ultracold atomic gases Also, the 
structures with energy in the second band and higher 
band-gaps warrant a deeper investigation, and it is con- 
ceivable that exact breather solutions could be found 
with structure similar to the ones that reappear in the 
dynamics. Beyond that, higher dimensional extensions 
would also be a challenging endeavor. 
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